load packages

library(tidyverse)
library(lme4)
library(lmerTest)
library(knitr)
library(MuMIn)
library(cowplot)

define color palettes

algorithm = c("#006989", "#FEC601", "#F43C13", "#00A5CF", "#00A878")
food = c("#A4C960", "#2A5C8C")
liking = c("#FF0000", "#F2AD00")
bid = c("#00A08A", "#F2AD00", "#F98400", "#FF0000")
dc_bw = readRDS("~/dc_bw.Rds")

load and tidy task data

task = read.csv("~/Documents/code/sanlab/DEV_scripts/fMRI/fx/multiconds/WTP/betaseries/events.csv", stringsAsFactors = FALSE) %>%
  mutate(bid = ifelse(bid == "NULL", NA, bid),
         bid = as.integer(bid),
         health = as.factor(health),
         liking = as.factor(liking)) %>%
  group_by(subjectID, wave, run) %>%
  mutate(trial = row_number()) %>%
  filter(wave == 1)

check responses

Check if wrong buttons were used (i.e., not 5-8)

  • DEV002 = wrong button box (i.e. 1-4) used
  • DEV007 = wrong button box (i.e. 1-4) used
  • DEV011 = code normally
  • DEV017 = exclude; can’t tell if they’re missed responses or incorrect placement of fingers
  • DEV019 = exclude; can’t tell if they’re missed responses or incorrect placement of fingers
  • DEV032 = code normally
  • DEV033 = incorrect placement of fingers; recode runs 1-3
  • DEV054 = exclude; technical error?
  • DEV061 = code normally
subs = task %>%
  mutate(bid = as.character(bid)) %>%
  group_by(subjectID, run, bid) %>%
  summarize(n = n()) %>%
  spread(bid, n) %>%
  mutate(messed = ifelse(!is.na(`2`), "yes", NA),
         messed = ifelse(is.na(`5`) & !is.na(`<NA>`), "yes", messed)) %>%
  filter(messed == "yes") %>% 
  ungroup() %>% 
  select(subjectID) %>% 
  unique()

task %>%
  mutate(bid = as.character(bid)) %>%
  group_by(subjectID, run, bid) %>%
  summarize(n = n()) %>%
  spread(bid, n) %>%
  mutate(messed = ifelse(!is.na(`2`), "yes", NA),
         messed = ifelse(is.na(`5`) & !is.na(`<NA>`), "yes", messed)) %>%
  filter(subjectID %in% subs$subjectID)

recode and exclude

Recoding
* DEV033: recode runs1-3, but if liking rating < 3, leave as missing

data.ex = task %>%
  mutate(bid = ifelse(subjectID == "DEV033" & !run == "run4", bid - 1, bid),
         bid = ifelse(subjectID == "DEV033" & !run == "run4" & is.na(bid) & liking_rating > 2, 8, bid),
         bid = (bid - 5) / 2) %>%
  group_by(subjectID, wave) %>%
  arrange(subjectID, run) %>%
  mutate(trial = row_number())

load striping info

striping = read.csv("~/Documents/code/sanlab/DEV_scripts/fMRI/betaseries/WTP/striping_QC.csv")

load mean intensity values

file_dir = "~/Documents/code/sanlab/DEV_scripts/fMRI/betaseries/WTP/dotProducts_WTP/"
file_pattern = "DEV[0-9]{3}_meanIntensity.txt"
file_list = list.files(file_dir, pattern = file_pattern)

intensities = data.frame()

for (file in file_list) {
  temp = tryCatch(read.table(file.path(file_dir,file), fill = TRUE) %>%
                    rename("subjectID" = V1,
                           "meanIntensity" = V3) %>%
                    extract(V2, "beta", "beta_([0-9]{4}).nii") %>%
                    mutate(beta = as.integer(beta)), error = function(e) message(file))
  intensities = rbind(intensities, temp)
  rm(temp)
}

load dot products

file_dir = "~/Documents/code/sanlab/DEV_scripts/fMRI/betaseries/WTP/dotProducts_WTP/"
file_pattern = "DEV[0-9]{3}_dotProducts.txt"
file_list = list.files(file_dir, pattern = file_pattern)

dots = data.frame()

for (file in file_list) {
  temp = tryCatch(read.table(file.path(file_dir,file), fill = TRUE) %>%
                    rename("subjectID" = V1,
                           "map" = V3,
                           "dotProduct" = V4) %>%
                    extract(V2, "beta", "beta_([0-9]{4}).nii") %>%
                    extract(map, "algorithm", "(.*)_.*.nii") %>%
                    mutate(beta = as.integer(beta)), error = function(e) message(file))
  dots = rbind(dots, temp)
  rm(temp)
}

join intensities and dots

  • recode trials with extreme intensities as NA
dots.merged = dots %>%
  left_join(., intensities, by = c("subjectID", "beta")) %>%
  group_by(subjectID, algorithm) %>%
  mutate(rownum = row_number())

# plot original
dots.merged %>%
  filter(algorithm == "logistic") %>%
  ggplot(aes(1, meanIntensity)) +
    geom_boxplot()

# assess extreme values and exclude when calculating SDs
dots.merged %>%
  filter(algorithm == "logistic") %>%
  arrange(meanIntensity)
dots.merged %>%
  filter(algorithm == "logistic") %>%
  arrange(-meanIntensity)
# recode outliers as NA
dots.merged = dots.merged %>%
  ungroup() %>%
  mutate(meanIntensity = ifelse(meanIntensity > 1 | meanIntensity < -1, NA, meanIntensity),
         median = median(meanIntensity, na.rm = TRUE),
         sd3 = 3*sd(meanIntensity, na.rm = TRUE),
         outlier = ifelse(meanIntensity > median + sd3 | meanIntensity < median - sd3, "yes", "no"),
         dotProduct = ifelse(outlier == "yes", NA, dotProduct))
  
# plot after
dots.merged %>%
  filter(algorithm == "logistic") %>%
  ggplot(aes(1, meanIntensity)) +
    geom_boxplot()

exclude outliers and standardize

  • standardize within algorithm and contrast
dots.ex = dots.merged %>%
  group_by(algorithm, subjectID) %>%
  mutate(trial = row_number()) %>%
  left_join(., striping, by = c("subjectID", "beta"))  %>%
  mutate(dotProduct = ifelse(!is.na(striping), NA, dotProduct)) %>%
  filter(!algorithm %in% c("ridge", "svm")) %>%
  group_by(subjectID, algorithm) %>% # standardize within sub and algorithm
  mutate(dotSTD = scale(dotProduct, center = FALSE)) 

merge data

Exclusions

  • Didn’t scan: DEV002, DEV007
  • MRI motion and data quality exclusions: DEV001, DEV020, DEV032, DEV047, DEV063, DEV067, DEV078
  • Button box exclusions: DEV017, DEV019, DEV054
  • Run exclusions: DEV028 (run1), DEV048 (run3), DEV064 (run1), DEV069 (run1)
data = left_join(dots.ex, data.ex, by = c("subjectID", "trial")) %>%
  filter(!subjectID %in% c("DEV002", "DEV007", "DEV001", "DEV020", "DEV032", "DEV047", "DEV063", "DEV067", "DEV078", "DEV017", "DEV019", "DEV054")) %>%
  filter(!(subjectID == "DEV028" & run == "run1") & !(subjectID == "DEV048" & run == "run3") & !(subjectID == "DEV064" & run == "run1") & !(subjectID == "DEV069" & run == "run1")) %>%
  ungroup() %>%
  ungroup() %>%
  mutate(algorithm = ifelse(algorithm == "reg_look", "regulate > look", 
                     ifelse(algorithm == "reg", "regulate > rest", algorithm)),
         liking = ifelse(liking_rating > 2, "liked",
                  ifelse(liking_rating < 3, "disliked", NA)))

descriptives

data %>%
  filter(!is.na(bid)) %>%
  ggplot(aes(bid, fill = liking)) +
    geom_histogram(position = "dodge") +
    facet_grid(~health) +
    scale_fill_manual(values = liking) +
    dc_bw

visualize

health

data %>%
  filter(!is.na(bid)) %>%
  ggplot(aes(algorithm, dotSTD, fill = health)) +
    stat_summary(fun.y = mean, geom = "bar", position = position_dodge(width = 0.95)) +
    stat_summary(fun.data = mean_cl_boot, geom = "errorbar", position = position_dodge(width = 0.95), width = 0) +
    scale_fill_manual(name = "", values = food) + 
    labs(y = "standarized regulation pattern expression value\n", x = "") + 
    dc_bw

data %>%
  filter(!is.na(bid)) %>%
  ggplot(aes(health, dotSTD)) +
    stat_summary(aes(group = subjectID), fun.y = mean, geom = "line", alpha = .1, size = .5) +
    stat_summary(aes(group = 1), fun.y = mean, geom = "line", size = .75) +
    stat_summary(aes(color = health), fun.data = mean_cl_boot, geom = "pointrange", width = 0, size = .75) + 
    facet_grid(~algorithm) +
    scale_color_manual(name = "", values = food) + 
    labs(y = "standarized regulation pattern expression value\n", x = "") + 
    dc_bw

health and liking

data %>%
  filter(!is.na(bid)) %>%
  ggplot(aes(algorithm, dotSTD, fill = health, alpha = liking)) +
    stat_summary(fun.y = mean, geom = "bar", position = position_dodge(width = 0.95)) +
    stat_summary(fun.data = mean_cl_boot, geom = "errorbar", position = position_dodge(width = 0.95), width = 0) +
    scale_fill_manual(name = "", values = food) + 
    scale_alpha_discrete(range = c(.6,1)) +
    labs(y = "standarized regulation pattern expression value\n", x = "") + 
    dc_bw

data %>%
  filter(!is.na(bid)) %>%
  ggplot(aes(liking, dotSTD)) +
    stat_summary(aes(group = interaction(subjectID, health), color = health), fun.y = mean, geom = "line", alpha = .1, size = .5) +
    stat_summary(aes(group = health, color = health), fun.y = mean, geom = "line", size = .75) +
    stat_summary(aes(color = health), fun.data = mean_cl_boot, geom = "pointrange", width = 0, size = .75) + 
    facet_grid(~algorithm) +
    scale_color_manual(name = "", values = food) + 
    coord_cartesian(ylim = c(-1.4, 1.2)) + 
    labs(y = "standarized regulation pattern expression value\n", x = "") + 
    dc_bw

bid value

data %>%
  filter(!is.na(bid)) %>%
  ggplot(aes(algorithm, dotSTD, fill = as.factor(bid))) +
    stat_summary(fun.y = mean, geom = "bar", position = position_dodge(width = 0.95)) +
    stat_summary(fun.data = mean_cl_boot, geom = "errorbar", position = position_dodge(width = 0.95), width = 0) +
    scale_fill_manual(name = "", values = bid) + 
    facet_grid(~health) + 
    labs(y = "standarized regulation pattern expression value\n", x = "") + 
    dc_bw

data %>%
  filter(!is.na(bid)) %>%
  group_by(bid, algorithm, health) %>%
  mutate(n.obs = n()) %>%
  ggplot(aes(as.factor(bid), dotSTD, color = algorithm)) +
    stat_summary(aes(group = algorithm), fun.y = mean, geom = "line") +
    stat_summary(fun.data = "mean_cl_boot", geom = "linerange") +
    stat_summary(aes(size = n.obs), fun.y = mean, geom = "point") +
    scale_color_manual(name = "", values = algorithm) + 
    scale_size(name = "", range = c(1,4)) + 
    facet_grid(~health) + 
    labs(x = "\nbid value", y = "standardized regulation pattern expression value\n") + 
    dc_bw

bid value and health

data %>%
  filter(!is.na(bid)) %>%
  group_by(bid, algorithm, health) %>%
  mutate(n.obs = n()) %>%
  ggplot(aes(as.factor(bid), dotSTD, color = health)) +
    stat_summary(aes(group = health), fun.y = mean, geom = "line") +
    stat_summary(fun.data = "mean_cl_boot", geom = "linerange") +
    stat_summary(aes(size = n.obs), fun.y = mean, geom = "point") +
    scale_color_manual(name = "", values = food) + 
    scale_size(name = "", range = c(1,4)) + 
    facet_grid(~algorithm) + 
    labs(x = "\nbid value", y = "standardized regulation pattern expression value\n") + 
    dc_bw

bid high v. low

all items

data %>%
  filter(!is.na(bid)) %>%
  mutate(bid.bin = ifelse(bid >= 1, "high", "low")) %>%
  ggplot(aes(bid.bin, dotSTD)) +
    stat_summary(aes(group = subjectID), fun.y = mean, geom = "line", alpha = .1, size = .5) +
    stat_summary(aes(group = 1), fun.y = mean, geom = "line", size = .75) +
    stat_summary(aes(color = bid.bin), fun.data = mean_cl_boot, geom = "pointrange", width = 0, size = .75) + 
    facet_grid(~algorithm) +
    scale_color_manual(name = "", values = liking) + 
    coord_cartesian(ylim = c(-1, 1.25)) + 
    labs(y = "standarized regulation pattern expression value\n", x = "") + 
    dc_bw

data %>%
  filter(!is.na(bid)) %>%
  mutate(bid.bin = ifelse(bid >= 1, "high", "low")) %>%
  ggplot(aes(bid.bin, dotSTD)) +
    stat_summary(aes(group = interaction(subjectID, health), color = health), fun.y = mean, geom = "line", alpha = .1, size = .5) +
    stat_summary(aes(group = health, color = health), fun.y = mean, geom = "line", size = .75) +
    stat_summary(aes(color = health), fun.data = mean_cl_boot, geom = "pointrange", width = 0, size = .75) + 
    facet_grid(~algorithm) +
    scale_color_manual(name = "", values = food) + 
    coord_cartesian(ylim = c(-1.25, 1.25)) + 
    labs(y = "standarized regulation pattern expression value\n", x = "") + 
    dc_bw

liked items only

data %>%
  filter(!is.na(bid)) %>%
  filter(liking == "liked") %>%
  mutate(bid.bin = ifelse(bid >= 1, "high", "low")) %>%
  ggplot(aes(bid.bin, dotSTD)) +
    stat_summary(aes(group = subjectID), fun.y = mean, geom = "line", alpha = .1, size = .5) +
    stat_summary(aes(group = 1), fun.y = mean, geom = "line", size = .75) +
    stat_summary(aes(color = bid.bin), fun.data = mean_cl_boot, geom = "pointrange", width = 0, size = .75) + 
    facet_grid(~algorithm) +
    scale_color_manual(name = "", values = liking) + 
    coord_cartesian(ylim = c(-1.5, 1.5)) + 
    labs(y = "standarized regulation pattern expression value\n", x = "") + 
    dc_bw

data %>%
  filter(!is.na(bid)) %>%
  filter(liking == "liked") %>%
  mutate(bid.bin = ifelse(bid >= 1, "high", "low")) %>%
  ggplot(aes(bid.bin, dotSTD)) +
    stat_summary(aes(group = interaction(subjectID, health), color = health), fun.y = mean, geom = "line", alpha = .1, size = .5) +
    stat_summary(aes(group = health, color = health), fun.y = mean, geom = "line", size = .75) +
    stat_summary(aes(color = health), fun.data = mean_cl_boot, geom = "pointrange", width = 0, size = .75) + 
    facet_grid(~algorithm) +
    scale_color_manual(name = "", values = food) + 
    coord_cartesian(ylim = c(-2, 2)) + 
    labs(y = "standarized regulation pattern expression value\n", x = "") + 
    dc_bw

RT

data %>%
  filter(!is.na(bid)) %>%
  ggplot(aes(rt, dotSTD, color = health)) +
    geom_smooth(method = "lm", alpha = .2) + 
    facet_grid(~algorithm) + 
    scale_color_manual(name = "", values = food) + 
    labs(x = "\nreaction time (seconds)", y = "standardized regulation pattern expression value\n") + 
    dc_bw

data %>%
  filter(!is.na(bid)) %>%
  ggplot(aes(rt, dotSTD, color = health, linetype = liking)) +
    geom_smooth(method = "lm", alpha = .1) + 
    facet_grid(~algorithm) + 
    scale_color_manual(name = "", values = food) + 
    labs(x = "\nreaction time (seconds)", y = "standardized regulation pattern expression value\n") + 
    dc_bw

data %>%
  filter(!is.na(bid)) %>%
  filter(algorithm == "logistic") %>%
  ggplot(aes(as.factor(bid), rt, color = health)) +
    stat_summary(aes(group = health), fun.y = mean, geom = "line") +
    stat_summary(fun.data = mean_cl_boot, geom = "pointrange", width = 0) +
    facet_grid(~liking) +
    scale_color_manual(values = food) +
    labs(x = "\nbid value", y = "reaction time (seconds)\n") + 
    dc_bw

data %>%
  filter(!is.na(bid)) %>%
  filter(algorithm == "logistic") %>%
  ggplot(aes(rt, dotSTD, color = as.factor(bid))) +
    geom_smooth(method = "lm", alpha = .1) + 
    facet_grid(~health) + 
    scale_color_manual(name = "bid", values = bid) + 
    labs(x = "\nreaction time (seconds)", y = "standardized regulation pattern expression value\n") + 
    dc_bw

data %>%
  filter(!is.na(bid)) %>%
  filter(algorithm == "logistic") %>%
  filter(liking == "liked") %>%
  ggplot(aes(rt, dotSTD, color = as.factor(bid))) +
    geom_smooth(method = "lm", alpha = .1) + 
    facet_grid(~health) + 
    scale_color_manual(name = "bid on liked foods", values = bid) + 
    labs(x = "\nreaction time (seconds)", y = "standardized regulation pattern expression value\n") + 
    dc_bw

individual diffs

all items

data %>%
  filter(!is.na(bid)) %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(meanPEV = mean(dotProduct, na.rm = TRUE)) %>%
  ggplot(aes(health, meanPEV, color = health)) + 
    geom_boxplot() +
    facet_wrap(~algorithm, scales = "free", ncol = 4) +
    scale_color_manual(values = food) +
    labs(y = "mean regulation pattern expression value\n") +
    dc_bw

data %>%
  filter(!is.na(bid)) %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
         meanBid = mean(bid, na.rm = TRUE)) %>%
  select(subjectID, algorithm, meanPEV, meanBid, health) %>%
  unique() %>%
  ggplot(aes(meanPEV, meanBid, color = health)) + 
    geom_point(alpha = .4) +
    geom_smooth(method = "lm", alpha = .2, fullrange = TRUE) +
    facet_wrap(~algorithm, scales = "free", ncol = 4) +
    scale_color_manual(values = food) +
    labs(x = "\nmean regulation pattern expression value", y = "mean bid value\n") +
    dc_bw

data %>%
  filter(!is.na(bid)) %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
         meanBid = mean(bid, na.rm = TRUE),
         health1 = health) %>%
  select(subjectID, algorithm, meanPEV, meanBid, health, health1) %>%
  unique() %>%
  spread(health, meanBid) %>%
  group_by(subjectID, algorithm) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(pref.health = healthy - unhealthy,
         pref.health.percent = ((healthy - unhealthy) / healthy) * 100) %>%
  spread(health1, meanPEV) %>%
  mutate(diff = unhealthy - healthy) %>%
  ggplot(aes(diff, pref.health.percent)) + 
    geom_point() +
    geom_smooth(method = "lm", se = .1) +
    facet_wrap(~algorithm, scales = "free", ncol = 4) +
    labs(x = "\ndifference in mean regulation pattern expression value (unhealthy - healthy)", y = "percent change in bid (healthy - uneahlty / healthy)") + 
    dc_bw

data %>%
  filter(!is.na(bid)) %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
         meanBid = mean(bid, na.rm = TRUE),
         health1 = health) %>%
  select(subjectID, algorithm, meanPEV, meanBid, health, health1) %>%
  unique() %>%
  spread(health, meanBid) %>%
  group_by(subjectID, algorithm) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(pref.health = healthy - unhealthy,
         pref.health.percent = ((healthy - unhealthy) / healthy) * 100) %>%
  spread(health1, meanPEV) %>%
  mutate(diff = unhealthy - healthy) %>% 
  ungroup() %>%
  nest(-algorithm) %>% 
  mutate(
    test = map(data, ~ cor.test(.x$pref.health.percent, .x$diff)),
    tidied = map(test, broom::tidy)
  ) %>% 
  unnest(tidied, .drop = TRUE)

liked items only

data %>%
  filter(!is.na(bid)) %>%
  filter(liking == "liked") %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(meanPEV = mean(dotProduct, na.rm = TRUE)) %>%
  ggplot(aes(health, meanPEV, color = health)) + 
    geom_boxplot() +
    facet_wrap(~algorithm, scales = "free", ncol = 4) +
    scale_color_manual(values = food) +
    labs(y = "mean regulation pattern expression value\n") +
    dc_bw

data %>%
  filter(!is.na(bid)) %>%
  filter(liking == "liked") %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
         meanBid = mean(bid, na.rm = TRUE)) %>%
  select(subjectID, algorithm, meanPEV, meanBid, health) %>%
  unique() %>%
  ggplot(aes(meanPEV, meanBid, color = health)) + 
    geom_point(alpha = .4) +
    geom_smooth(method = "lm", alpha = .2, fullrange = TRUE) +
    facet_wrap(~algorithm, scales = "free", ncol = 4) +
    scale_color_manual(values = food) +
    labs(x = "\nmean regulation pattern expression value", y = "mean bid value\n") +
    dc_bw

data %>%
  filter(!is.na(bid)) %>%
  filter(liking == "liked") %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
         meanBid = mean(bid, na.rm = TRUE),
         health1 = health) %>%
  select(subjectID, algorithm, meanPEV, meanBid, health, health1) %>%
  unique() %>%
  spread(health, meanBid) %>%
  group_by(subjectID, algorithm) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(pref.health = healthy - unhealthy,
         pref.health.percent = ((healthy - unhealthy) / healthy) * 100) %>%
  spread(health1, meanPEV) %>%
  mutate(diff = unhealthy - healthy) %>%
  ggplot(aes(diff, pref.health.percent)) + 
    geom_point() +
    geom_smooth(method = "lm", se = .1) +
    facet_wrap(~algorithm, scales = "free", ncol = 4) +
    labs(x = "\ndifference in mean regulation pattern expression value (unhealthy - healthy)", y = "percent change in bid (healthy - uneahlty / healthy)") + 
    dc_bw

data %>%
  filter(!is.na(bid)) %>%
  filter(liking == "liked") %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
         meanBid = mean(bid, na.rm = TRUE),
         health1 = health) %>%
  select(subjectID, algorithm, meanPEV, meanBid, health, health1) %>%
  unique() %>%
  spread(health, meanBid) %>%
  group_by(subjectID, algorithm) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(pref.health = healthy - unhealthy,
         pref.health.percent = ((healthy - unhealthy) / healthy) * 100) %>%
  spread(health1, meanPEV) %>%
  mutate(diff = unhealthy - healthy) %>% 
  ungroup() %>%
  nest(-algorithm) %>% 
  mutate(
    test = map(data, ~ cor.test(.x$pref.health.percent, .x$diff)),
    tidied = map(test, broom::tidy)
  ) %>% 
  unnest(tidied, .drop = TRUE)

percent

all items

data %>%
  filter(!is.na(bid)) %>%
  filter(!subjectID == "DEV037") %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
         meanBid = mean(bid, na.rm = TRUE),
         regulation = ifelse(dotSTD > 0, 1, 0),
         sumReg = sum(regulation, na.rm = TRUE),
         n = n(), 
         percentReg = (sumReg / n) * 100) %>%
  select(subjectID, algorithm, meanPEV, meanBid, percentReg, health) %>%
  unique() %>%
  ggplot(aes(percentReg, meanBid, color = health)) + 
    geom_point(alpha = .4) +
    geom_smooth(method = "lm", alpha = .2, fullrange = TRUE) +
    facet_wrap(~algorithm, scales = "free", ncol = 4) +
    scale_color_manual(values = food) +
    labs(x = "\npercent regulation", y = "mean bid value\n") +
    dc_bw

data %>%
  filter(!is.na(bid)) %>%
  filter(!subjectID == "DEV037") %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
         meanBid = mean(bid, na.rm = TRUE),
         regulation = ifelse(dotSTD > 0, 1, 0),
         sumReg = sum(regulation, na.rm = TRUE),
         n = n(), 
         percentReg = (sumReg / n) * 100) %>%
  select(subjectID, algorithm, meanPEV, meanBid, percentReg, health) %>%
  unique() %>%
  ggplot(aes(percentReg, meanPEV, color = health)) + 
    geom_point(alpha = .4) +
    geom_smooth(method = "lm", alpha = .2, fullrange = TRUE) +
    facet_wrap(~algorithm, scales = "free", ncol = 4) +
    scale_color_manual(values = food) +
    labs(x = "\npercent regulation", y = "mean regulation pattern expression valu\n") +
    dc_bw

data %>%
  filter(!is.na(bid)) %>%
  filter(!subjectID == "DEV037") %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
         meanBid = mean(bid, na.rm = TRUE),
         regulation = ifelse(dotSTD > 0, 1, 0),
         sumReg = sum(regulation, na.rm = TRUE),
         n = n(), 
         percentReg = (sumReg / n) * 100) %>%
  select(subjectID, algorithm, meanPEV, meanBid, percentReg, health) %>%
  unique() %>%
  group_by(algorithm, health) %>%
  mutate(medReg = median(percentReg),
         binReg = ifelse(percentReg >= median(percentReg), "high", "low")) %>%
  ggplot(aes(meanPEV, meanBid, color = health, linetype = binReg)) + 
    geom_point(alpha = .2) +
    geom_smooth(method = "lm", alpha = .2, fullrange = TRUE) +
    facet_wrap(~algorithm, scales = "free", ncol = 4) +
    scale_color_manual(values = food) +
    labs(x = "\nmean regulation pattern expression value", y = "mean bid value\n") +
    dc_bw

data %>%
  filter(!is.na(bid)) %>%
  filter(!subjectID == "DEV037") %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
         meanBid = mean(bid, na.rm = TRUE),
         regulation = ifelse(dotSTD > 0, 1, 0),
         sumReg = sum(regulation, na.rm = TRUE),
         n = n(), 
         percentReg = (sumReg / n) * 100,
         health1 = health,
         health2 = health) %>%
  select(subjectID, algorithm, meanPEV, meanBid, percentReg, health, health1, health2) %>%
  unique() %>%
  spread(health, meanBid) %>%
  group_by(subjectID, algorithm) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(pref.health = healthy - unhealthy,
         pref.health.percent = ((healthy - unhealthy) / healthy) * 100) %>%
  spread(health1, meanPEV) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(diff = unhealthy - healthy) %>%
  spread(health2, percentReg) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(diffReg = unhealthy - healthy) %>%
  unique() %>%
  ggplot(aes(diffReg, pref.health.percent)) + 
    geom_point() +
    geom_smooth(method = "lm", se = .1) +
    facet_wrap(~algorithm, scales = "free", ncol = 4) +
    labs(x = "\ndifference percent regulation (unhealthy - healthy)", y = "percent change in bid (healthy - uneahlty / healthy)") + 
    dc_bw

data %>%
  filter(!is.na(bid)) %>%
  filter(!subjectID == "DEV037") %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
         meanBid = mean(bid, na.rm = TRUE),
         regulation = ifelse(dotSTD > 0, 1, 0),
         sumReg = sum(regulation, na.rm = TRUE),
         n = n(), 
         percentReg = (sumReg / n) * 100,
         health1 = health,
         health2 = health) %>%
  select(subjectID, algorithm, meanPEV, meanBid, percentReg, health, health1, health2) %>%
  unique() %>%
  spread(health, meanBid) %>%
  group_by(subjectID, algorithm) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(pref.health = healthy - unhealthy,
         pref.health.percent = ((healthy - unhealthy) / healthy) * 100) %>%
  spread(health1, meanPEV) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(diff = unhealthy - healthy) %>%
  spread(health2, percentReg) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(diffReg = unhealthy - healthy) %>%
  ungroup() %>%
  nest(-algorithm) %>% 
  mutate(
    test = map(data, ~ cor.test(.x$pref.health.percent, .x$diffReg)),
    tidied = map(test, broom::tidy)
  ) %>% 
  unnest(tidied, .drop = TRUE)
data %>%
  filter(!is.na(bid)) %>%
  filter(!subjectID == "DEV037") %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
         meanBid = mean(bid, na.rm = TRUE),
         regulation = ifelse(dotSTD > 0, 1, 0),
         sumReg = sum(regulation, na.rm = TRUE),
         n = n(), 
         percentReg = (sumReg / n) * 100,
         health1 = health,
         health2 = health) %>%
  select(subjectID, algorithm, meanPEV, meanBid, percentReg, health, health1, health2) %>%
  unique() %>%
  spread(health, meanBid) %>%
  group_by(subjectID, algorithm) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(pref.health = healthy - unhealthy,
         pref.health.percent = ((healthy - unhealthy) / healthy) * 100) %>%
  spread(health1, meanPEV) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(diff = unhealthy - healthy) %>%
  spread(health2, percentReg) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(diffReg = unhealthy - healthy) %>%
  unique() %>%
  ggplot(aes(diffReg, diff)) + 
    geom_point() +
    geom_smooth(method = "lm", se = .1) +
    facet_wrap(~algorithm, scales = "free", ncol = 4) +
    labs(x = "\ndifference percent regulation (unhealthy - healthy)", y = "difference in mean regulation pattern expression value (unhealthy - healthy") + 
    dc_bw

data %>%
  filter(!is.na(bid)) %>%
  filter(!subjectID == "DEV037") %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
         meanBid = mean(bid, na.rm = TRUE),
         regulation = ifelse(dotSTD > 0, 1, 0),
         sumReg = sum(regulation, na.rm = TRUE),
         n = n(), 
         percentReg = (sumReg / n) * 100,
         health1 = health,
         health2 = health) %>%
  select(subjectID, algorithm, meanPEV, meanBid, percentReg, health, health1, health2) %>%
  unique() %>%
  spread(health, meanBid) %>%
  group_by(subjectID, algorithm) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(pref.health = healthy - unhealthy,
         pref.health.percent = ((healthy - unhealthy) / healthy) * 100) %>%
  spread(health1, meanPEV) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(diff = unhealthy - healthy) %>%
  spread(health2, percentReg) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(diffReg = unhealthy - healthy) %>%
  ungroup() %>%
  nest(-algorithm) %>% 
  mutate(
    test = map(data, ~ cor.test(.x$diff, .x$diffReg)),
    tidied = map(test, broom::tidy)
  ) %>% 
  unnest(tidied, .drop = TRUE)

liked items only

data %>%
  filter(!is.na(bid)) %>%
  filter(liking == "liked") %>%
  filter(!subjectID == "DEV037") %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
         meanBid = mean(bid, na.rm = TRUE),
         regulation = ifelse(dotSTD > 0, 1, 0),
         sumReg = sum(regulation, na.rm = TRUE),
         n = n(), 
         percentReg = (sumReg / n) * 100) %>%
  select(subjectID, algorithm, meanPEV, meanBid, percentReg, health) %>%
  unique() %>%
  ggplot(aes(percentReg, meanBid, color = health)) + 
    geom_point(alpha = .4) +
    geom_smooth(method = "lm", alpha = .2, fullrange = TRUE) +
    facet_wrap(~algorithm, scales = "free", ncol = 4) +
    scale_color_manual(values = food) +
    labs(x = "\npercent regulation", y = "mean bid value\n") +
    dc_bw

data %>%
  filter(!is.na(bid)) %>%
  filter(liking == "liked") %>%
  filter(!subjectID == "DEV037") %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
         meanBid = mean(bid, na.rm = TRUE),
         regulation = ifelse(dotSTD > 0, 1, 0),
         sumReg = sum(regulation, na.rm = TRUE),
         n = n(), 
         percentReg = (sumReg / n) * 100) %>%
  select(subjectID, algorithm, meanPEV, meanBid, percentReg, health) %>%
  unique() %>%
  ggplot(aes(percentReg, meanPEV, color = health)) + 
    geom_point(alpha = .4) +
    geom_smooth(method = "lm", alpha = .2, fullrange = TRUE) +
    facet_wrap(~algorithm, scales = "free", ncol = 4) +
    scale_color_manual(values = food) +
    labs(x = "\npercent regulation", y = "mean regulation pattern expression valu\n") +
    dc_bw

data %>%
  filter(!is.na(bid)) %>%
  filter(liking == "liked") %>%
  filter(!subjectID == "DEV037") %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
         meanBid = mean(bid, na.rm = TRUE),
         regulation = ifelse(dotSTD > 0, 1, 0),
         sumReg = sum(regulation, na.rm = TRUE),
         n = n(), 
         percentReg = (sumReg / n) * 100) %>%
  select(subjectID, algorithm, meanPEV, meanBid, percentReg, health) %>%
  unique() %>%
  group_by(algorithm, health) %>%
  mutate(medReg = median(percentReg),
         binReg = ifelse(percentReg >= median(percentReg), "high", "low")) %>%
  ggplot(aes(meanPEV, meanBid, color = health, linetype = binReg)) + 
    geom_point(alpha = .2) +
    geom_smooth(method = "lm", alpha = .2, fullrange = TRUE) +
    facet_wrap(~algorithm, scales = "free", ncol = 4) +
    scale_color_manual(values = food) +
    labs(x = "\nmean standardized regulation pattern expression value", y = "mean bid value\n") +
    dc_bw

data %>%
  filter(!is.na(bid)) %>%
  filter(liking == "liked") %>%
  filter(!subjectID == "DEV037") %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
         meanBid = mean(bid, na.rm = TRUE),
         regulation = ifelse(dotSTD > 0, 1, 0),
         sumReg = sum(regulation, na.rm = TRUE),
         n = n(), 
         percentReg = (sumReg / n) * 100,
         health1 = health,
         health2 = health) %>%
  select(subjectID, algorithm, meanPEV, meanBid, percentReg, health, health1, health2) %>%
  unique() %>%
  spread(health, meanBid) %>%
  group_by(subjectID, algorithm) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(pref.health = healthy - unhealthy,
         pref.health.percent = ((healthy - unhealthy) / healthy) * 100) %>%
  spread(health1, meanPEV) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(diff = unhealthy - healthy) %>%
  spread(health2, percentReg) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(diffReg = unhealthy - healthy) %>%
  unique() %>%
  ggplot(aes(diffReg, pref.health.percent)) + 
    geom_point() +
    geom_smooth(method = "lm", se = .1) +
    facet_wrap(~algorithm, scales = "free", ncol = 4) +
    labs(x = "\ndifference percent regulation (unhealthy - healthy)", y = "percent change in bid (healthy - uneahlty / healthy)") + 
    dc_bw

data %>%
  filter(!is.na(bid)) %>%
  filter(liking == "liked") %>%
  filter(!subjectID == "DEV037") %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
         meanBid = mean(bid, na.rm = TRUE),
         regulation = ifelse(dotSTD > 0, 1, 0),
         sumReg = sum(regulation, na.rm = TRUE),
         n = n(), 
         percentReg = (sumReg / n) * 100,
         health1 = health,
         health2 = health) %>%
  select(subjectID, algorithm, meanPEV, meanBid, percentReg, health, health1, health2) %>%
  unique() %>%
  spread(health, meanBid) %>%
  group_by(subjectID, algorithm) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(pref.health = healthy - unhealthy,
         pref.health.percent = ((healthy - unhealthy) / healthy) * 100) %>%
  spread(health1, meanPEV) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(diff = unhealthy - healthy) %>%
  spread(health2, percentReg) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(diffReg = unhealthy - healthy) %>%
  ungroup() %>%
  nest(-algorithm) %>% 
  mutate(
    test = map(data, ~ cor.test(.x$pref.health.percent, .x$diffReg)),
    tidied = map(test, broom::tidy)
  ) %>% 
  unnest(tidied, .drop = TRUE)
data %>%
  filter(!is.na(bid)) %>%
  filter(liking == "liked") %>%
  filter(!subjectID == "DEV037") %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
         meanBid = mean(bid, na.rm = TRUE),
         regulation = ifelse(dotSTD > 0, 1, 0),
         sumReg = sum(regulation, na.rm = TRUE),
         n = n(), 
         percentReg = (sumReg / n) * 100,
         health1 = health,
         health2 = health) %>%
  select(subjectID, algorithm, meanPEV, meanBid, percentReg, health, health1, health2) %>%
  unique() %>%
  spread(health, meanBid) %>%
  group_by(subjectID, algorithm) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(pref.health = healthy - unhealthy,
         pref.health.percent = ((healthy - unhealthy) / healthy) * 100) %>%
  spread(health1, meanPEV) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(diff = unhealthy - healthy) %>%
  spread(health2, percentReg) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(diffReg = unhealthy - healthy) %>%
  unique() %>%
  ggplot(aes(diffReg, diff)) + 
    geom_point() +
    geom_smooth(method = "lm", se = .1) +
    facet_wrap(~algorithm, scales = "free", ncol = 4) +
    labs(x = "\ndifference percent regulation (unhealthy - healthy)", y = "difference in mean regulation pattern expression value (unhealthy - healthy") + 
    dc_bw

data %>%
  filter(!is.na(bid)) %>%
  filter(liking == "liked") %>%
  filter(!subjectID == "DEV037") %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
         meanBid = mean(bid, na.rm = TRUE),
         regulation = ifelse(dotSTD > 0, 1, 0),
         sumReg = sum(regulation, na.rm = TRUE),
         n = n(), 
         percentReg = (sumReg / n) * 100,
         health1 = health,
         health2 = health) %>%
  select(subjectID, algorithm, meanPEV, meanBid, percentReg, health, health1, health2) %>%
  unique() %>%
  spread(health, meanBid) %>%
  group_by(subjectID, algorithm) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(pref.health = healthy - unhealthy,
         pref.health.percent = ((healthy - unhealthy) / healthy) * 100) %>%
  spread(health1, meanPEV) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(diff = unhealthy - healthy) %>%
  spread(health2, percentReg) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(diffReg = unhealthy - healthy) %>%
  ungroup() %>%
  nest(-algorithm) %>% 
  mutate(
    test = map(data, ~ cor.test(.x$diff, .x$diffReg)),
    tidied = map(test, broom::tidy)
  ) %>% 
  unnest(tidied, .drop = TRUE)

restraint

# load qualtrics
cred_file_location = '~/credentials.yaml.DEFAULT'
sid_column_name = 'Login|qid'
survey_name_filter = 'DEV Session 1 Surveys'
sid_pattern = 'DEV[0-9]{3}'
exclude_sid = '99|DEV737' # subject IDs to exclude

# load credential file
credentials = scorequaltrics::creds_from_file(cred_file_location)

# filter
surveysAvail = scorequaltrics::get_surveys(credentials)
surveysFiltered = filter(surveysAvail, grepl(survey_name_filter, SurveyName))

# get data
surveys_long = scorequaltrics::get_survey_data(surveysFiltered,
                                               credentials, 
                                               pid_col = sid_column_name) %>%
  mutate(Login = gsub("Dev", "DEV", Login),
         Login = gsub("dev", "DEV", Login),
         Login = ifelse(grepl("^[0-9]{3}$", Login), paste0("DEV", Login), Login),
         Login = ifelse(Login == "DEVO55", "DEV055", Login))

# check for repeats
repeat.subs = surveys_long %>%
  filter(item == "Finished" & value == "1") %>%
  filter(!grepl(exclude_sid, Login) & !is.na(Login)) %>%
  group_by(survey_name, Login) %>%
  summarize(n = n()) %>%
  filter(n > 1) %>%
  spread(survey_name, n)

# surveys_long %>%
#   filter(item == "StartDate") %>%
#   filter(Login %in% repeat.subs$Login) %>%
#   group_by(survey_name, Login) %>%
#   mutate(n = n()) %>%
#   filter(n > 1)

# filter out repeats
surveys = surveys_long %>%
  filter(Login %in% unique(data$subjectID)) %>%
  filter(!qid %in% c("R_1M01CRpEgQ9Sjzx", "R_3JfnLZ2XhekmvvF", "R_RUXzgKp7Sne7865")) %>%
  select(-qid) %>%
  rename("subjectID" = Login)

# get and score restraint scale
restraint = surveys %>%
  filter(grepl("RS", item)) %>%
  mutate(value = ifelse(value == "", NA, value),
         value = as.integer(value)) %>%
  group_by(subjectID) %>%
  summarize(restraint = sum(value, na.rm = TRUE))

all items

data %>%
  filter(!is.na(bid)) %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(mean.sub = mean(dotProduct, na.rm = TRUE),
         mean.bid = mean(bid, na.rm = TRUE),
         health1 = health) %>%
  select(subjectID, algorithm, mean.sub, mean.bid, health, health1) %>%
  unique() %>%
  spread(health, mean.bid) %>%
  group_by(subjectID, algorithm) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(pref.health = healthy - unhealthy,
         pref.health.percent = ((healthy - unhealthy) / healthy) * 100) %>%
  spread(health1, mean.sub) %>%
  mutate(diff = unhealthy - healthy) %>%
  left_join(., restraint) %>%
  ggplot(aes(diff, restraint)) + 
    geom_point() +
    geom_smooth(method = "lm", se = .1) +
    facet_wrap(~algorithm, scales = "free", ncol = 3) +
    labs(x = "\ndifference in mean regulation pattern expression value (unhealthy - healthy)", y = "restraint score\n") +
    dc_bw

data %>%
  filter(!is.na(bid)) %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(mean.sub = mean(dotProduct, na.rm = TRUE),
         mean.bid = mean(bid, na.rm = TRUE),
         health1 = health) %>%
  select(subjectID, algorithm, mean.sub, mean.bid, health, health1) %>%
  unique() %>%
  spread(health, mean.bid) %>%
  group_by(subjectID, algorithm) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(pref.health = healthy - unhealthy,
         pref.health.percent = ((healthy - unhealthy) / healthy) * 100) %>%
  spread(health1, mean.sub) %>%
  mutate(diff = unhealthy - healthy) %>%
  left_join(., restraint) %>% 
  ungroup() %>%
  nest(-algorithm) %>% 
  mutate(
    test = map(data, ~ cor.test(.x$restraint, .x$diff)),
    tidied = map(test, broom::tidy)
  ) %>% 
  unnest(tidied, .drop = TRUE)
data %>%
  filter(!is.na(bid)) %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
         meanBid = mean(bid, na.rm = TRUE),
         regulation = ifelse(dotSTD > 0, 1, 0),
         sumReg = sum(regulation, na.rm = TRUE),
         n = n(), 
         percentReg = (sumReg / n) * 100,
         health1 = health,
         health2 = health) %>%
  select(subjectID, algorithm, meanPEV, meanBid, percentReg, health, health1, health2) %>%
  unique() %>%
  spread(health, meanBid) %>%
  group_by(subjectID, algorithm) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(pref.health = healthy - unhealthy,
         pref.health.percent = ((healthy - unhealthy) / healthy) * 100) %>%
  spread(health1, meanPEV) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(diff = unhealthy - healthy) %>%
  spread(health2, percentReg) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(diffReg = unhealthy - healthy) %>%
  left_join(., restraint) %>%
  ggplot(aes(diffReg, restraint)) + 
    geom_point() +
    geom_smooth(method = "lm", se = .1) +
    facet_wrap(~algorithm, scales = "free", ncol = 3) +
    labs(x = "\ndifference in percent regulation (unhealthy - healthy)", y = "restraint score\n") +
    dc_bw

data %>%
  filter(!is.na(bid)) %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
         meanBid = mean(bid, na.rm = TRUE),
         regulation = ifelse(dotSTD > 0, 1, 0),
         sumReg = sum(regulation, na.rm = TRUE),
         n = n(), 
         percentReg = (sumReg / n) * 100,
         health1 = health,
         health2 = health) %>%
  select(subjectID, algorithm, meanPEV, meanBid, percentReg, health, health1, health2) %>%
  unique() %>%
  spread(health, meanBid) %>%
  group_by(subjectID, algorithm) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(pref.health = healthy - unhealthy,
         pref.health.percent = ((healthy - unhealthy) / healthy) * 100) %>%
  spread(health1, meanPEV) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(diff = unhealthy - healthy) %>%
  spread(health2, percentReg) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(diffReg = unhealthy - healthy) %>%
  left_join(., restraint) %>% 
  ungroup() %>%
  nest(-algorithm) %>% 
  mutate(
    test = map(data, ~ cor.test(.x$restraint, .x$diffReg)),
    tidied = map(test, broom::tidy)
  ) %>% 
  unnest(tidied, .drop = TRUE)

liked items only

data %>%
  filter(!is.na(bid)) %>%
  filter(liking == "liked") %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(mean.sub = mean(dotProduct, na.rm = TRUE),
         mean.bid = mean(bid, na.rm = TRUE),
         health1 = health) %>%
  select(subjectID, algorithm, mean.sub, mean.bid, health, health1) %>%
  unique() %>%
  spread(health, mean.bid) %>%
  group_by(subjectID, algorithm) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(pref.health = healthy - unhealthy,
         pref.health.percent = ((healthy - unhealthy) / healthy) * 100) %>%
  spread(health1, mean.sub) %>%
  mutate(diff = unhealthy - healthy) %>%
  left_join(., restraint) %>%
  ggplot(aes(diff, restraint)) + 
    geom_point() +
    geom_smooth(method = "lm", se = .1) +
    facet_wrap(~algorithm, scales = "free", ncol = 3) +
    labs(x = "\ndifference in mean regulation pattern expression value (regulate - look)", y = "restraint score\n") +
    dc_bw

data %>%
  filter(!is.na(bid)) %>%
  filter(liking == "liked") %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(mean.sub = mean(dotProduct, na.rm = TRUE),
         mean.bid = mean(bid, na.rm = TRUE),
         health1 = health) %>%
  select(subjectID, algorithm, mean.sub, mean.bid, health, health1) %>%
  unique() %>%
  spread(health, mean.bid) %>%
  group_by(subjectID, algorithm) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(pref.health = healthy - unhealthy,
         pref.health.percent = ((healthy - unhealthy) / healthy) * 100) %>%
  spread(health1, mean.sub) %>%
  mutate(diff = unhealthy - healthy) %>%
  left_join(., restraint) %>% 
  ungroup() %>%
  nest(-algorithm) %>% 
  mutate(
    test = map(data, ~ cor.test(.x$restraint, .x$diff)),
    tidied = map(test, broom::tidy)
  ) %>% 
  unnest(tidied, .drop = TRUE)
data %>%
  filter(!is.na(bid)) %>%
  filter(liking == "liked") %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
         meanBid = mean(bid, na.rm = TRUE),
         regulation = ifelse(dotSTD > 0, 1, 0),
         sumReg = sum(regulation, na.rm = TRUE),
         n = n(), 
         percentReg = (sumReg / n) * 100,
         health1 = health,
         health2 = health) %>%
  select(subjectID, algorithm, meanPEV, meanBid, percentReg, health, health1, health2) %>%
  unique() %>%
  spread(health, meanBid) %>%
  group_by(subjectID, algorithm) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(pref.health = healthy - unhealthy,
         pref.health.percent = ((healthy - unhealthy) / healthy) * 100) %>%
  spread(health1, meanPEV) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(diff = unhealthy - healthy) %>%
  spread(health2, percentReg) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(diffReg = unhealthy - healthy) %>%
  left_join(., restraint) %>%
  ggplot(aes(diffReg, restraint)) + 
    geom_point() +
    geom_smooth(method = "lm", se = .1) +
    facet_wrap(~algorithm, scales = "free", ncol = 3) +
    labs(x = "\ndifference in percent regulation (unhealthy - healthy)", y = "restraint score\n") +
    dc_bw

data %>%
  filter(!is.na(bid)) %>%
  filter(liking == "liked") %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
         meanBid = mean(bid, na.rm = TRUE),
         regulation = ifelse(dotSTD > 0, 1, 0),
         sumReg = sum(regulation, na.rm = TRUE),
         n = n(), 
         percentReg = (sumReg / n) * 100,
         health1 = health,
         health2 = health) %>%
  select(subjectID, algorithm, meanPEV, meanBid, percentReg, health, health1, health2) %>%
  unique() %>%
  spread(health, meanBid) %>%
  group_by(subjectID, algorithm) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(pref.health = healthy - unhealthy,
         pref.health.percent = ((healthy - unhealthy) / healthy) * 100) %>%
  spread(health1, meanPEV) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(diff = unhealthy - healthy) %>%
  spread(health2, percentReg) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  mutate(diffReg = unhealthy - healthy) %>%
  left_join(., restraint) %>% 
  ungroup() %>%
  nest(-algorithm) %>% 
  mutate(
    test = map(data, ~ cor.test(.x$restraint, .x$diffReg)),
    tidied = map(test, broom::tidy)
  ) %>% 
  unnest(tidied, .drop = TRUE)

correlations

all items

cors = data %>%
  filter(!is.na(bid)) %>%
  group_by(subjectID, algorithm, health, liking) %>%
  mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
         meanBid = mean(bid, na.rm = TRUE)) %>%
  gather(variable, value, meanPEV, meanBid) %>%
  select(subjectID, health, liking, algorithm, variable, value) %>%
  unique() %>%
  unite("health", c("liking", "variable", "health")) %>%
  ungroup() %>%
  select(subjectID, algorithm, health, value) %>%
  group_by(algorithm) %>%
  do({
    health.spread = spread(., health, value)
    cors = cor(health.spread[,-c(1:2)], use = "pairwise.complete.obs") %>%
      as.data.frame() %>%
      mutate(algorithm = health.spread$algorithm[[1]],
             health = colnames(health.spread)[-c(1:2)])
  })


cors %>%
  reshape2::melt() %>%
  ggplot(aes(health, variable, fill = value)) +
    geom_tile(color = "white") +
    scale_fill_gradientn(name = "correlation\n", colors = c("#3B9AB2", "white", "#F21A00"), limits = c(-1, 1), breaks = c(-1, 0, 1)) + 
    geom_text(aes(label = round(value, 2)), size = 3) + #family = "Futura Medium"
    #geom_text(aes(label = significant), size = 4, family = "Futura Medium", nudge_x = .3, nudge_y = .2) + 
    facet_wrap(~algorithm) +
    labs(x = "", y = "") + 
    dc_bw +
    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))

cors = data %>%
  filter(!is.na(bid)) %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
         meanBid = mean(bid, na.rm = TRUE)) %>%
  gather(variable, value, meanPEV, meanBid) %>%
  select(subjectID, health, algorithm, variable, value) %>%
  unique() %>%
  unite("health", c("variable", "health")) %>%
  ungroup() %>%
  select(subjectID, algorithm, health, value) %>%
  group_by(algorithm) %>%
  do({
    health.spread = spread(., health, value)
    cors = cor(health.spread[,-c(1:2)], use = "pairwise.complete.obs") %>%
      as.data.frame() %>%
      mutate(algorithm = health.spread$algorithm[[1]],
             health = colnames(health.spread)[-c(1:2)])
  })


cors %>%
  reshape2::melt() %>%
  ggplot(aes(health, variable, fill = value)) +
    geom_tile(color = "white") +
    scale_fill_gradientn(name = "correlation\n", colors = c("#3B9AB2", "white", "#F21A00"), limits = c(-1, 1), breaks = c(-1, 0, 1)) + 
    geom_text(aes(label = round(value, 2)), size = 3) + #family = "Futura Medium"
    #geom_text(aes(label = significant), size = 4, family = "Futura Medium", nudge_x = .3, nudge_y = .2) + 
    facet_wrap(~algorithm) +
    labs(x = "", y = "") + 
    dc_bw +
    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))

liked items only

cors = data %>%
  filter(!is.na(bid)) %>%
  filter(liking == "liked") %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
         meanBid = mean(bid, na.rm = TRUE)) %>%
  gather(variable, value, meanPEV, meanBid) %>%
  select(subjectID, health, algorithm, variable, value) %>%
  unique() %>%
  unite("health", c("variable", "health")) %>%
  ungroup() %>%
  select(subjectID, algorithm, health, value) %>%
  group_by(algorithm) %>%
  do({
    health.spread = spread(., health, value)
    cors = cor(health.spread[,-c(1:2)], use = "pairwise.complete.obs") %>%
      as.data.frame() %>%
      mutate(algorithm = health.spread$algorithm[[1]],
             health = colnames(health.spread)[-c(1:2)])
  })


cors %>%
  reshape2::melt() %>%
  ggplot(aes(health, variable, fill = value)) +
    geom_tile(color = "white") +
    scale_fill_gradientn(name = "correlation\n", colors = c("#3B9AB2", "white", "#F21A00"), limits = c(-1, 1), breaks = c(-1, 0, 1)) + 
    geom_text(aes(label = round(value, 2)), size = 3) + #family = "Futura Medium"
    #geom_text(aes(label = significant), size = 4, family = "Futura Medium", nudge_x = .3, nudge_y = .2) + 
    facet_wrap(~algorithm) +
    labs(x = "", y = "") + 
    dc_bw +
    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))

unhealthy items only

cors = data %>%
  filter(!is.na(bid)) %>%
  filter(health == "unhealthy") %>%
  group_by(subjectID, algorithm) %>%
  mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
         meanBid = mean(bid, na.rm = TRUE)) %>%
  gather(variable, value, meanPEV, meanBid) %>%
  select(subjectID, algorithm, variable, value) %>%
  unique() %>%
  ungroup() %>%
  select(subjectID, algorithm, variable, value) %>%
  group_by(algorithm) %>%
  do({
    health.spread = spread(., variable, value)
    cors = cor(health.spread[,-c(1:2)], use = "pairwise.complete.obs", method = "pearson") %>%
      as.data.frame() %>%
      mutate(algorithm = health.spread$algorithm[[1]],
             health = colnames(health.spread)[-c(1:2)])
  })


cors %>%
  reshape2::melt() %>%
  ggplot(aes(health, variable, fill = value)) +
    geom_tile(color = "white") +
    scale_fill_gradientn(name = "correlation\n", colors = c("#3B9AB2", "white", "#F21A00"), limits = c(-1, 1), breaks = c(-1, 0, 1)) + 
    geom_text(aes(label = round(value, 2)), size = 3) + #family = "Futura Medium"
    #geom_text(aes(label = significant), size = 4, family = "Futura Medium", nudge_x = .3, nudge_y = .2) + 
    facet_wrap(~algorithm) +
    labs(x = "", y = "") + 
    dc_bw +
    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))

liked unhealthy items only

cors = data %>%
  filter(!is.na(bid)) %>%
  filter(liking == "liked" & health == "unhealthy") %>%
  group_by(subjectID, algorithm) %>%
  mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
         meanBid = mean(bid, na.rm = TRUE)) %>%
  gather(variable, value, meanPEV, meanBid) %>%
  select(subjectID, algorithm, variable, value) %>%
  unique() %>%
  ungroup() %>%
  select(subjectID, algorithm, variable, value) %>%
  group_by(algorithm) %>%
  do({
    health.spread = spread(., variable, value)
    cors = cor(health.spread[,-c(1:2)], use = "pairwise.complete.obs", method = "pearson") %>%
      as.data.frame() %>%
      mutate(algorithm = health.spread$algorithm[[1]],
             health = colnames(health.spread)[-c(1:2)])
  })


cors %>%
  reshape2::melt() %>%
  ggplot(aes(health, variable, fill = value)) +
    geom_tile(color = "white") +
    scale_fill_gradientn(name = "correlation\n", colors = c("#3B9AB2", "white", "#F21A00"), limits = c(-1, 1), breaks = c(-1, 0, 1)) + 
    geom_text(aes(label = round(value, 2)), size = 3) + #family = "Futura Medium"
    #geom_text(aes(label = significant), size = 4, family = "Futura Medium", nudge_x = .3, nudge_y = .2) + 
    facet_wrap(~algorithm) +
    labs(x = "", y = "") + 
    dc_bw +
    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))

data %>%
  filter(!is.na(health)) %>%
  filter(health == "unhealthy" & liking == "liked") %>%
  group_by(subjectID, algorithm, health) %>%
  mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
         meanBid = mean(bid, na.rm = TRUE)) %>%
  select(subjectID, algorithm, meanPEV, meanBid) %>%
  unique() %>% 
  ungroup() %>%
  nest(-algorithm) %>% 
  mutate(
    test = map(data, ~ cor.test(.x$meanBid, .x$meanPEV)),
    tidied = map(test, broom::tidy)
  ) %>% 
  unnest(tidied, .drop = TRUE)

specification curve analysis

health

models

# set na.action for dredge
options(na.action = "na.fail")

# tidy data
data.sca = data %>%
  select(subjectID, trial, algorithm, dotSTD, bid, health) %>%
  spread(algorithm, dotSTD) %>%
  na.omit()

# run full model
models = lmer(bid ~ health*logistic + health*`regulate > look` + health*`regulate > rest` + (1 | subjectID), data = data.sca)

# run all possible nested models
models.sca = dredge(models, rank = "AIC", extra = "BIC") %>%
  select(AIC, delta, BIC, df, logLik, weight, `(Intercept)`, health, everything())

# set AIC for null model you want to compare model AIC values to
null = models.sca %>%
  arrange(df) %>%
  filter(health == "+" & df == 4)

plot

# tidy for plotting
plot.data = models.sca %>%
  arrange(AIC) %>%
  mutate(specification = row_number(),
         better.fit = ifelse(AIC == null$AIC, "equal",
                      ifelse(AIC < null$AIC, "yes", "no")))

order = plot.data %>%
  arrange(AIC) %>%
  mutate(better.fit.num = ifelse(better.fit == "yes", 1, 0)) %>%
  gather(variable, value, -c(AIC, delta, BIC, df, logLik, weight, better.fit.num, specification, better.fit)) %>%
  filter(!is.na(value)) %>%
  group_by(variable) %>%
  mutate(order = sum(better.fit.num)) %>%
  select(variable, order) %>%
  unique()

# variables included in model
variable.names = names(select(plot.data, -starts_with("better"), -specification, -AIC, -BIC, -df, -logLik, -delta, -weight, -`(Intercept)`))

# plot top panel
a = plot.data %>%
  ggplot(aes(specification, AIC, color = better.fit)) +
    geom_point(shape = "|", size = 4, alpha = .75) +
    geom_hline(yintercept = null$AIC, linetype = "dashed", color = "#5BBCD6") +
    scale_color_manual(values = c("#5BBCD6", "black", "#F43C13")) +
    labs(x = "", y = "AIC\n") +
    theme_minimal(base_size = 14) +
    theme(legend.title = element_text(size = 1),
          legend.text = element_text(size = 9),
          axis.text = element_text(color = "black"),
          axis.line = element_line(colour = "black"),
          legend.position = "none",
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          panel.background = element_blank())

# set plotting order for variables based on number of times it's included in better fitting models
b = plot.data %>%
  gather(variable, value, eval(variable.names)) %>%
  left_join(., order, by = "variable") %>%
  mutate(value = ifelse(!is.na(value), "|", ""),
         variable = gsub("`regulate > look`", "regulate > look", variable),
         variable = gsub("`regulate > rest`", "regulate > rest", variable),
         variable = gsub("regulate > look:health", "health:regulate > look", variable),
         variable = gsub("regulate > rest:health", "health:regulate > rest", variable),
         variable = gsub("health:", "health  x  ", variable)) %>%
  ggplot(aes(specification, reorder(variable, order), color = better.fit)) +
    geom_text(aes(label = value), alpha = .75) +
    scale_color_manual(values = c("#5BBCD6", "black", "#F43C13")) +
    labs(x = "\nspecification number", y = "variables\n") +
    theme_minimal(base_size = 14) +
    theme(legend.title = element_text(size = 1),
          legend.text = element_text(size = 9),
          axis.text = element_text(color = "black"),
          axis.line = element_line(colour = "black"),
          legend.position = "none",
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          panel.background = element_blank())

(plot = plot_grid(a, b, ncol = 1, align = "v"))

liking and health

models

# set na.action for dredge
options(na.action = "na.fail")

# tidy data
data.sca = data %>%
  select(subjectID, trial, algorithm, dotSTD, bid, health, liking) %>%
  spread(algorithm, dotSTD) %>%
  na.omit()

# run full model
models = lmer(bid ~ liking*logistic + liking*`regulate > look` + liking*`regulate > rest` + health*logistic + health*`regulate > look` + health*`regulate > rest` +(1 | subjectID), data = data.sca)

# run all possible nested models
models.sca = dredge(models, rank = "AIC", extra = "BIC") %>%
  select(AIC, delta, BIC, df, logLik, weight, `(Intercept)`, liking, everything())

# set AIC for null model you want to compare model AIC values to
null = models.sca %>%
  arrange(df) %>%
  filter(liking == "+" & df == 4)

plot

# tidy for plotting
plot.data = models.sca %>%
  arrange(AIC) %>%
  mutate(specification = row_number(),
         better.fit = ifelse(AIC == null$AIC, "equal",
                      ifelse(AIC < null$AIC, "yes", "no")))

order = plot.data %>%
  arrange(AIC) %>%
  mutate(better.fit.num = ifelse(better.fit == "yes", 1, 0)) %>%
  gather(variable, value, -c(AIC, delta, BIC, df, logLik, weight, better.fit.num, specification, better.fit)) %>%
  filter(!is.na(value)) %>%
  group_by(variable) %>%
  mutate(order = sum(better.fit.num)) %>%
  select(variable, order) %>%
  unique()

# variables included in model
variable.names = names(select(plot.data, -starts_with("better"), -specification, -AIC, -BIC, -df, -logLik, -delta, -weight, -`(Intercept)`))

# plot top panel
a = plot.data %>%
  ggplot(aes(specification, AIC, color = better.fit)) +
    geom_point(shape = "|", size = 4, alpha = .75) +
    geom_hline(yintercept = null$AIC, linetype = "dashed", color = "#5BBCD6") +
    scale_color_manual(values = c("#5BBCD6", "black", "#F43C13")) +
    labs(x = "", y = "AIC\n") +
    theme_minimal(base_size = 14) +
    theme(legend.title = element_text(size = 1),
          legend.text = element_text(size = 9),
          axis.text = element_text(color = "black"),
          axis.line = element_line(colour = "black"),
          legend.position = "none",
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          panel.background = element_blank())

# set plotting order for variables based on number of times it's included in better fitting models
b = plot.data %>%
  gather(variable, value, eval(variable.names)) %>%
  left_join(., order, by = "variable") %>%
  mutate(value = ifelse(!is.na(value), "|", ""),
         variable = gsub("`regulate > look`", "regulate > look", variable),
         variable = gsub("`regulate > rest`", "regulate > rest", variable),
         variable = gsub("regulate > look:liking", "liking:regulate > look", variable),
         variable = gsub("regulate > rest:liking", "liking:regulate > rest", variable),
         variable = gsub("liking:", "liking  x  ", variable),
         variable = gsub("regulate > look:health", "health:regulate > look", variable),
         variable = gsub("regulate > rest:health", "health:regulate > rest", variable),
         variable = gsub("health:", "health  x  ", variable)) %>%
  ggplot(aes(specification, reorder(variable, order), color = better.fit)) +
    geom_text(aes(label = value), alpha = .75) +
    scale_color_manual(values = c("#5BBCD6", "black", "#F43C13")) +
    labs(x = "\nspecification number", y = "variables\n") +
    theme_minimal(base_size = 14) +
    theme(legend.title = element_text(size = 1),
          legend.text = element_text(size = 9),
          axis.text = element_text(color = "black"),
          axis.line = element_line(colour = "black"),
          legend.position = "none",
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          panel.background = element_blank())

(plot = plot_grid(a, b, ncol = 1, align = "v"))